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Abstract. This paper presents a theoretical and numerical investigation of the following prac- 
tical question: Should the time difference quotients of the snapshots be used to generate the proper 
orthogonal decomposition basis functions? The answer to this question is important, since some 
published numerical studies use the time difference quotients, whereas other numerical studies do 
not. The criterion used in this paper to answer this question is the rate of convergence of the error 
of the reduced order model with respect to the number of proper orthogonal decomposition basis 
functions. Two cases are considered: the noJDQ case, in which the snapshot difference quotients are 
not used, and the DQ case, in which the snapshot difference quotients are used. The error estimates 
suggest that the convergence rates in the C°(L 2 ) -norm and in the C^-ff 1 )-norm are optimal for 
the DQ case, but suboptimal for the no_DQ case. The convergence rates in the L 2 (// 1 )-norm are 
optimal for both the DQ case and the no_DQ case. Numerical tests are conducted on the heat 
equation and on the Burgers equation. The numerical results support the conclusions drawn from 
the theoretical error estimates. Overall, the theoretical and numerical results strongly suggest that, 
in order to achieve optimal pointwise in time rates of convergence with respect to the number of 
proper orthogonal decomposition basis functions, one should use the snapshot difference quotients. 
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1. Introduction. This paper addresses the following question: "Should the time 
difference quotients (DQs) of the snapshots be used in the generation of the Proper 
Orthogonal Decomposition (POD) basis functions?" 

We emphasize that this is an important question. There are two schools of 
thought: one uses the DQs [13j El E| , the other does not [17 ] [18 ] 14] [21 ] . 

To our knowledge, the first instance in which the snapshot DQs were incorporated 
in the generation of the POD basis functions was the pioneering paper of Kunisch and 
Volkwcin [13] . In that report, the authors considered two types of errors for a general 
parabolic equation: the time discretization errors and the POD truncation errors. 
They argued that one needs to include the temporal difference quotients in the set of 
snapshots; otherwise, the error will be suboptimal, containing an extra factor (see 
Remark 1 in |13j). Thus, the motivation for using the temporal difference quotients 
was purely theoretical. In numerical investigations, however, the authors reported 
contradictory findings: in |13) . the use of the DQs did not improve the quality of the 
reduced-order model; in [9], however, they did. Kunisch and Volkwein used again the 
snapshot DQs when they considered the Navier-Stokes equations [14] . 

The snapshot DQs were also used in the Discrete Empirical Interpolation Method 
(DEIM) of Chaturantabut and Sorensen (6] [5]. Their motivation, however, was dif- 
ferent from that of Kunisch and Volkwein. Indeed, the authors considered in [5] [5] 
a general, nonlinear system of equations of the form y' = f(y, t). In the set of 
snapshots, they included not only the state variables y, but also the nonlinear snap- 
shots f(y,t). They further noted (see page 48 in [5]) that, since f(y,£) = y' and 
(y™ +1 — y")/Ai ~ y', this is similar to including the temporal difference quotients, 
as done in [l3l [14] . 
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To our knowledge, the first reports in which the DQs were not used were [17j for 
the heat equation and [32] for the Navier-Stokes equations. Chapelle et al. g] used 
a different approach that did not use the DQs either. This approach employed the 
L 2 projection instead of the standard H 1 projection used in, e.g., [13l[T4j- Further 
improvements to the approach used in [4] (as well as that used in p~3l [14] ) were made 
by Singler in [21"] . 

From the above discussion, it is clear that the question whether the snapshot DQs 
should be included or not in the set of snapshots is important. To our knowledge, this 
question is still open. This report represents a first step in answering this question. 
All our discussion will be centered around the heat equation, although most (if not all) 
of it could be extended to general parabolic equations in a straightforward manner. 

From a theoretical point of view, the only motivation for using the snapshot DQs 
was given in Remark 1 in [13] . The main point of this remark is the following: In the 
error analysis of the evolution equation, to approximate Ut(t n ), the time derivative of 
the exact solution u evaluated at time t n , one usually uses the continuous difference 
quotient du(t n ) :— — — 27"^ • ^ approximate the continuous difference quotient 
du(t n ) in the POD space, one naturally uses the POD difference quotient du r {t n ) := 

— 27"^* " • We assume that u r (t n+1 ) is an optimal approximation for u(t n+1 ) 
and that u r (t n ) is an optimal approximation for u(t n ), where the optimality is with 
respect to r and At. Then, it would appear that, with respect with At, du r (t n ) is 
a suboptimal approximation for du{t n ), because of the At in the denominator of the 
two difference quotients. 

Although the argument above, used in Remark 1 in [13] to motivate the inclusion 
of the snapshot DQs in the derivation of the POD basis, seems natural, we point 
out that this issue should be treated more carefully. Indeed, in the finite element 
approximation of parabolic equations, it is well known that the DQs duh(t n ) := 

— 2^"" fel ' t - are actually optimal (with respect with At) approximations of the 
continuous difference quotients du(t n ) (see, e.g., [HI [20]) . Thus, since the POD 
and finite element approximations are similar (both use a Galerkin projection in the 
spatial discretization) , one could question the validity of the argument used in Remark 
1. We emphasize that we are not claiming that the above argument is not valid; we 
are merely pointing out that a rigorous numerical analysis is needed before drawing 
any conclusions. 

Our preliminary numerical studies indicate that not using the DQs docs not yield 
suboptimal (with respect to At) error estimates. For the heat equation (see Section^] 
for details regarding the parameters used) , we monitor the rates of convergence with 
respect to At for the POD reduced-order model. We consider two cases: when the DQs 
are used in the generation of the POD basis (the corresponding results are denoted 
by DQ), and when the DQs are not used in the generation of the POD basis (the 
corresponding results are denoted by no-DQ). The errors are listed in Table ITTT1 and 
plotted in Figure H~T1 with associated linear regressions (LR). Both no-DQ and DQ 
approaches obtain optimal approximation order 0(At 2 ) in the L 2 norm. 

The rest of the paper is organized as follows: In Section[2l we sketch the derivation 
of the POD reduced order model. In Section^ we carefully derive the error estimates 
for the POD reduced order model. We focus on the rates of convergence with respect 
to r, the number of POD basis functions. In Section^ we present numerical results for 
two test problems: the heat equation and the Burgers equation. Finally, in Section [5j 
we draw several conclusions. 



DIFFERENCE QUOTIENTS IN PROPER ORTHOGONAL DECOMPOSITION 3 




At 



Fig. 1.1: Heat equation. Plots of the errors in the L 2 norm with respect to the time 
step At when the DQs are used (denoted by DQ) and when the DQs are not used 
(denoted by no-DQ). 



Table 1.1: Errors of the no-DQ and DQ approach when At varies. 



At 


no-DQ 


DQ 


r 


£ L 2 




r 


£ L 2 




2.00e-01 


6 


3.71e-02 


9.26c-01 


6 


3.71e-02 


9.26e-01 


1.00e-01 


11 


1.27e-02 


5.81c-01 


11 


1.27e-02 


5.81e-01 


5.00e-02 


21 


2.99e-03 


1.97e-01 


21 


2.99e-03 


1.97e-01 


2.50e-02 


41 


6.53e-04 


3.81e-02 


41 


6.53e-04 


3.81e-02 


1.00e-02 


59 


1.03e-04 


1.15e-02 


88 


1.03e-04 


1.15e-02 



2. Proper Orthogonal Decomposition Reduced Order Modeling. In this 
section, we sketch the derivation of the standard POD Galcrkin reduced-order model 
for the heat equation. For a detailed presentation of this reduced-order model, in- 
cluding more general settings, the reader is referred to [8l H~2l [22] . 

Let X := Hq(CI), where ft is the computational domain. Let u(-,t) e X, t £ [0,T] 
be the weak solution of the weak formulation of the heat equation: 

(ut,«) + i/(Vu,V«) = (/,«) MveX. (2.1) 

Given the time instances t±,. . .,tN £ [0, T], we consider the following two ensembles 
of snapshots: 

yno.DQ := span { u (. , to ) f . . . f u (. , tjv )} , (2.2) 

V DQ :=span{u(;t ),...,u(;t N ),du(;t 1 ),...,du(;t N )} , (2.3) 

where du{t n ) := ; n = 1, . . . , N are the time difference quotients (DQs). 

The two ensembles of snapshots correspond to the two cases investigated in this paper: 
(i) with the DQs not included in the snapshots (i.e., V no - D< ^); and (ii) with the DQs 
included in the snapshots (i.e., V D< ^). To simplify the presentation, we denote by V 
both sets of snapshots and use the specific notation (i.e., y n °- D Q or V D ®) only when 
is necessary. Let dim V = d. The POD seeks a low-dimensional basis {(pi, . . . , ^r}, 
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with r <^ d, which optimally approximates the input collection. Specifically, the POD 
basis satisfies 



1 N 



3=1 



(2.4) 



L 2 



subject to the conditions that (ipi,ipj)i,2 — 5y, 1 < i,j < r. 
Definition 2.1. TTie ierm 

r 

V inter P(x,t) :=u{x,t)-Y,H^t),vAx)) L . *>i(*) 

3=1 



(2.5) 



wiH 6e denoted as the POD interpolation error. 

Remark 2.1. Note that the H 1 norm can also be used to generate the POD basis 
JMWM- In this case, ||r/ m * er P|| L 2 ~ || V?7 m * erp || L 2 . Thus, the two cases considered in 
this paper (i.e., V = V no - DQ orV = V DQ ) yield error estimates that have the same 
convergence rates with respect to r. 

In order to solve (|2.4j) . we consider the eigenvalue problem 



K v = A v . 



(2.6) 



where K <S R NxN , with = — (u(-,tj),u(-,ti)) L2 , is the snapshot correlation 

matrix, Ai > A2 > . . . > > are the positive eigenvalues, and Vk, k = 1, . . . , d are 
the associated eigenvectors. It can then be shown (see, e.g., [HI HI]), that the solution 
of (|2.4p is given by 



1 N 



(2.7) 



where (vk)j is the j-th component of the eigenvector Vk- It can also be shown |13) 
that the following POD approximation property holds: 



1 N 

-Y 

j=i 



3=1 



j=r+l 



— Y 

IN + 1 f-^ 



3=1 



L- 



3=1 



L' 2 



if V = y"°- D o , (2.8) 
(2.9) 



^ A, HV = V DQ . 

j=r+l 



The approximation property (|2.8|) - ()2.9|) represents the relationship between the av- 
erage of the square of the L 2 -norm of the interpolation error -A. ^2 i=1 \\rf nterp {-, 
and the sum of the eigenvalues of the POD modes that are not included in the POD 
reduced order model. 

In what follows, we will use the notation X r — spanj^i, if2, ■ ■ ■ , ip r } ■ To derive 
the POD reduced-order model for the heat equation (|2.ip . we employ the Galerkin 
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truncation, which yields the following approximation u r £ X r of u: 

r 

u r {x,t):=Y,aj{t)vAx)- (2-10) 

3=1 

Plugging (|2.10[) into (|2.1[) and multiplying by test functions in X r C X yields the 
POD Galerkin reduced order model (POD-G-ROM): 

(u rt t,V r ) + v(Vu r ,Vv r ) = (f,V r ) Vl> r € X T . (2.11) 

The main advantage of the POD-G-ROM (|2.1ip over a straightforward finite element 
discretization of (|2.ip is clear - the computational cost of the former is dramatically 
lower than that of the latter. 

3. Error Estimates. In this section, we prove estimates for the error ||u — u r \\ , 
where u is the solution of the weak formulation of the heat equation (|2.1[) and u r is 
the solution of the POD-G-ROM (|2.1ip . In our theoretical analysis, we consider two 
cases, depending on the type of snapshots used in the derivation of the POD basis: 
Case I: V = V DQ (i.e., with the DQs); and Case II: V = V no - DQ (i.e., without the 
DQs). The main goal of this paper is to investigate whether Case I, Case II, or 
both Case I and Case II, yield error estimates that are optimal with respect to r. 
The optimality with respect to r is given by the following error estimates: 

||«-«r|| <C\\r, interp \\, (3.1) 
\\X7(u-u r )\\ <C\\Vrf nterp \\ , (3.2) 

where rf nterp is the POD interpolation error defined in (|2 . 5[) . 

We emphasize that ||77 m * erp || and || Vrf nterp \\ scale differently with respect to r: 
The scaling of ||7y mter P|| is given by the approximation property (|2.8I) - (|2.9[) . (In fact, 
only the average of ||7y mte, "P|| is presented. Since we will exclusively consider the 
continuous in time formulation, without loss of generality, we assume that that the 
pointwisc in time version of (|2.8j) - (|2.9j) is also valid.) The scaling of ||V?7 m * erp || is 
not given by the approximation property (|2.8|) - (|2.9j) . To derive such an estimate, we 
use the fact that the interpolation error lives in a finite dimensional space, i.e., the 
space spanned by the snapshots. Using an inverse estimate similar to that presented 
in Lemma 13.11 but for the entire space of snapshots (of dimension d) , we get the 
following estimate: 

||VT? <nte ^|| £a < C mv {d) \\ V mterp \\ L , , (3.3) 

where Ci nv (d) is the constant in the inverse estimate in Lemma 13.11 Following the 
discussion in Remark I3.R we conclude that the scaling of |j Vn mterp \\ is of lower order 
with respect to r than the scaling of ||77" l * erp || . Thus, if the error analysis yields 
estimates of the form 

||«-«r|| <C||Vrr* erp ||, (3.4) 

these estimates will be called suboptimal with respect to r. 

In this section, we investigate the optimality question from a theoretical point 
of view, by monitoring the dependency of the error estimates on r. In Section 21 we 
investigate the same question from a numerical point of view. 
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Wc note that we perform the error analysis only for the semidiscretization of the 
POD-G-ROM (|2.1ip . In fact, in this semidiscretization, wc only consider the error 
component corresponding to the POD truncation. Of course, in practical numerical 
simulations, the semidiscretization also has a spatial component (e.g., due to the 
finite clement discretization - see Section [4j. Furthermore, when considering the 
full discretization, the error also has a time discretization component (e.g., due to 
the time stepping algorithm - see Section |4]) . All these error components should be 
included in a rigorous error analysis of the discretization of the POD-G-ROM (|2.11[) 
(see, e.g., |17[[TUl[TTj ). For clarity of presentation, however, we only consider the error 
component corresponding to the POD truncation. In what follows, we will show that 
this is sufficient for answering the question asked in the title. 

We start by introducing some notation and we list several results that will be used 
throughout this section. For clarity, wc will denote by C a generic constant that can 
depend on all the parameters in the system, except on N (the number of snapshots), 
d (the dimension of the set of snapshots, V), and r (the number of POD modes used 
in the POD-G-ROMfl2X[])). 

We note that, since the POD basis is computed in the L 2 -norm (see (|2.4j) ). the 
POD mass matrix M r £ W xr with Mij = (<Pj,y>i) is the identity matrix. Thus, the 
POD inverse estimate that was proven in [TU] becomes: 

Lemma 3.1 (POD Inverse Estimate). Let S r <G R rxr with = (V<pj, V<fi) be 
the POD stiffness matrix and \\ ■ U2 denote the matrix 2-norm. Then, for all v r £ X r , 
the following POD inverse estimate holds: 

||Vv r || ia < C inv (r) \\v r \\ L 2 , (3.5) 

where C inv (r) := v /||S , r || 2 . 

Remark 3.1 (POD Inverse Estimate Scalings). Since the r dependency of the 
error estimates presented in this section will be carefully monitored, we try to get some 
insight into the scalings of the constant Ci nv (r) in (|3.5p . i-c, the scalings of 1 1 jSV- 1 1 2 
with respect to r, the dimension of the POD basis. We note that, since the POD 
basis significantly varies from test case to test case, it would be difficult to derive a 
general scaling of\\S r \\2- We emphasize, however, that when the underlying system is 
homogeneous (i.e., invariant to spatial translations), the POD basis is identical to the 
Fourier basis (see, e.g., JE/)- I n this case, it is easy to derive the scalings of\\S r \\2 with 
respect to r. Without loss of generality, we assume that the computational domain is 
[0, 1] C K 1 and that the boundary conditions are homogeneous Dirichlet. In this case, 
the Fourier basis functions are given by the following formula: <Pj(x) = sin(j7rx). It 
is a simple calculation to show that the matrix S r is diagonal and that its diagonal 
entries are given by the following formula: 

Sjj = J (j tt) 2 sin 2 (j 7T x) dx = - (j tt) 2 . (3.6) 

It is easy to see that, in the n-dimensional case, formula (|3.6p becomes 

% = /(i7r) 2 " sm 2 (jirx)dx= \{jv? n - (3.7) 

Since the POD stiffness matrix S r is symmetric, its matrix 2-norm is given by ||SV||2 = 
A mol , where A mQX is the largest eigenvalue of S r . Thus, in the n-dimensional case, 
we have 

\\S r \\ 2 = ^(r-K) 2n = 0(r 2n ). (3.8) 
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Thus, we conclude that, when the underlying system is homogeneous, the 2-norm of 
the POD stiffness matrix S r scales as 0(r 2n ), where n is the spatial dimension. 

As mentioned at the beginning of the remark, for general (non-homogeneous) sys- 
tems is would be hard to derive theoretical scalings. The numerical tests in Section [7J 
however, seem to confirm the theoretical scaling in (|3.8[) . 

For the heat equation and the Burgers equation (see Section^for details regarding 
the parameters used), we monitor the scaling of \\S r \\ with respect to r for the POD 
reduced order model. We consider two cases: when the DQs are used in the generation 
of the POD basis (the corresponding results are denoted by DQ), and when the DQs 
are not used in the generation of the POD basis (the corresponding results are denoted 
by no-DQ). The scalings are plotted in Fiaure [?Tl\ It is seen that both no-DQ and DQ 
approaches yields scalings that are similar to the theoretical scaling (|3.8I) predicted for 
the homogeneous flow fields (i.e., when the POD basis reduces to the Fourier basis). 
The only exception seems to be for the Burgers equation in the DQ case. In all cases, 
however, the scaling \\S r \\2 = 0(r a ), where a is a positive constant, seems to be valid. 



no-DQ 

- LR order 2.25 
DQ 

- LR order 1 .95 



° no-DQ 
LR order 2.89 

□ DQ 
LR order 0.67 



Fig. 3.1: Heat equation (left); Burgers equation (right). Plots of the scalings of \\S r \\ 
with respect to r when the DQs are used (denoted by DQ) and when the DQs are not 
used (denoted by no-DQ). 



After these preliminaries, we are ready to derive the error estimates. The error 
analysis will proceed along the same lines as the error analysis of the finite element 
discretization [16]. The main difference between the finite element and the POD 
settings is that the finite element approximation property is global I16j , whereas 
the POD approximation property is local, i.e., it is only valid around the snapshots 
(see (|2.8j) - (|2.9[) 'l. This difference between the finite element and POD settings has ma- 
jor implications on the error analysis. This implications will be emphasized through- 
out this section. 

We start by considering the error equation: 

(e t ,w r ) + i/(Ve,Vw r ) =0 Vu r G X r , (3.9) 
where e := u — u r is the error. The error is split into two parts: 

e = u — u r = (u — w r ) — (u r — w r ) = r\ — 4> r , (3.10) 

where w r is an arbitrary function in X r , r\ := u — w r , and <j) r := u r — w r . Using this 
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decomposition in the error equation (|3.9[) . we get 

(<j>r,t,v r ) + v (V(j) r ,Vv r ) = {r] t ,v r ) + v{V-q, Vv r ). (3-11) 

The analysis proceeds by using (|3.11j) to show that 

Hr\\<C\\vl (3.12) 

where C is a generic constant that does not depend on r. Using the triangle inequality, 
one then gets 

||e|| < NI + lhM < (1 + C)||r?||. (3.13) 

Since w r was chosen arbitrarily and since the LHS does not depend on w r , we can 
take the infimum over w r in (|3.13j) and get the following error estimate: 

II 6 II < (1 + C) inf \\u-w r \\. (3.14) 

w r £X T 

At this stage, one invokes the approximability property of the space X r in (|3.14|l and 
concludes that the error estimate p,14[) is optimal with respect to r. That is, the error 
in (|3.14p is, up to a constant, the interpolation error in X r . 

There are a lot of variations of the error analysis sketched above, but most of the 
existing error analyses follow the above path. The main variations are the following: 
(i) the choice of the arbitrary function w r £ X r ; (ii) the norm || • |J used above; and 
(iii) the choice of the test function used in the error equation to get (|3.12p . 

In the remainder of this section, we investigate whether error estimates that are 
optimal with respect to r can be obtained with or without including the DQs in the 
set of snapshots. To this end, in Section [37X1 we consider the case in which the DQs 
are included in the set of snapshots (i.e., V = V D ®). Then, in Section [5721 we consider 
the case in which the DQs are not included in the set of snapshots (i.e., V = y no - D Qy 

3.1. Case I (V = V D ®). The standard approach used to prove error estimates 
in this case is to use the Ritz projection JT3J [14j [TUJ [IT] . We note that this is the 
standard approach used in the finite element context [16] . As pointed out in [23] , 
the Ritz projection was first used by Wheeler in |24j to obtain optimal error estimates 
for the finite element discretization of parabolic problems. 

We start by choosing w r := R r {u) in (|3.10|) . where R r {u) is the Ritz projection of 
u, given by: 

(V(tt - Rr(u)), Vu r ) = y Vr eX r . (3.15) 
Using (I3.15p . (|3.1ip becomes 



{<f>r,t, v r ) + v (V0 r , Viv) = (r}t,v r ) + v(^j^^hp). (3.16) 

It is the cancelation of the last term on the RHS of (|3.16p that yields optimal error 
estimates. We let v r := 4> r in (|3 . 1 6|) . and then we apply the Cauchy-Schwarz inequality 
to the remaining term on the RHS: 

^ll^-|| 2 + ^||V0 r || 2 <j|ry t j|||0 r ||. (3.17) 
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We rewrite the first term on the LHS of (|3.17p as 

lj t Ur\\ 2 = Ur\\^Url (3.18) 

Wc apply the Poincare-Friedrichs inequality to the second term on the LHS of (|3.17[) : 

H|V0 r || 2 > CV||<A r || 2 . (3.19) 

We note that the Poincare-Friedrichs inequality C||?;|| 2 < ||Vu|| 2 holds for every func- 
tion v in the continuous space Hq(CI), and, in particular, for <j> r £ X r C X = iJg(Q) 
(see equation (3) in [IS])- Thus, the constant C in (|3.19[) does not depend on r. Using 
(|3~Tg)) and (j3~TOl) in ([3~T?1) . we get 

lj t Hr\\+CpUr\\<\\Vtl (3.20) 
Using Gronwall's lemma in p.20p . we get 

||<M*)II <e- 2C ^|l0r(O)||+2 f e-^^^h^Wds. (3.21) 

Jo 

Using f|3 . 1 3[) . the first term on the RHS of (|3 . 2 1 1) can be estimated as follows: 

||^(0)|| < ||e(0)|| + h(0)||. (3.22) 

Thus, (|3~2T|) becomes 

\\Mt)\\ < e~ 2Cvt (||e(0)|| + ||„(0)||) + 2 £ e -^(*-) || % ( S )|| da. (3.23) 
Applying the triangle inequality, just as in (|3.13[) . we get 

||e(t)|| < Ht)\\ +er 2C » t f||e(0)|| + h(0)||") + 2 f ' e~ 2Gv{t - s) \\vt(s)\\ds. (3.24) 



Estimate (|3. 24[) shows that, as long as the Ritz projection (|3.15[) yields estimates 
for ||?7(t)|| (including for t — 0) and ||?7t(t)|| that are optimal with respect to r, the 
estimates for the POD error e are also optimal with respect to r. 

Remark 3.2. In the finite element context, both \\r] Rltz (t)\\ and ||i7t i,t *(i)|| are 
optimal (with respect to h, the mesh size). In the POD context, however, this is not 
that clear. To the best of the authors' knowledge, the state-of-the-art regarding the 
Ritz projection in a POD context is given in the pioneering paper of Kunisch and 
Volkwein jlSj. Since the Ritz projection plays such an important role in this paper, 
we summarize below the results in 

We start by emphasizing that the full discretization (i.e., in time and in the POD 
space) is used in whereas we only use the semidiscretization (i.e., in the POD 

space). Therefore, a significant difference between the two settings is that i] £ X r in 

(since only the discrete time instances tk are considered), whereas r) G X in our 
setting. 

The main result in regarding the Ritz projection is Lemma 3, which, in our 
notation, states the following: 

d 

liVry^l! 2 < !|S r j| 2 J2 V ( 3 - 25 ) 

j=r+l 
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For clarity of presentation, in (|3.25p . we have not included the constants that do not 
depend on r. Thus, the following relationship between the Ritz projection error and 
the POD interpolation error holds: 

||V77 ffit l <C^/\S4 2 \\rf nter v\\ . (3.26) 

The scaling in (|3.26p suggests that the Ritz projection yields optimal error estimates 
with respect to r in the H 1 -seminorm (see Q3.3[) ). 

We emphasize that Lemma 3 in U3\j does not include any bounds for the L 2 norm 
of n Rltz , i.e., for ||i7 *||. This is in clear contrast with the finite element context, in 
which \\r) Rttz \\ is estimated by the usual duality argument (the Aubin-Nitsche "trick"). 
Using a duality argument, however, is challenging in the POD context, since any 
auxiliary dual problem would not necessarily inherit the POD approximation prop- 
erty (|2.8[) or (|2.9[) . To the best of the authors ' knowledge, such a duality argument 
has never been used in a POD context. 

We emphasize that not being able to use a duality argument in the Ritz projection 
to get error estimates that are optimal with respect to r has significant consequences 
in the error analysis. Indeed, in the proof of Theorem 7 in J 13}/ (the error estimate 
for the backward Euler time discretization), to estimate the \\r) R,ltz \\ error component 
in equations (27a) and (27b), the authors use the WWij^^W estimate given in Lemma 
3 and the Poincare-Friedrichs inequality given in equation (3). Since the Poincare- 
Friedrichs constant does not depend on r, we conclude that \\f] mtz \\ and IIVt^* 2 !! 
have the same order. This, in turn, suggests that \\r] Rltz \\ is suboptimal with respect 
to r. We note that the same approach (i.e., Lemma 4 and the Poincare-Friedrichs 
inequality) is used in f!3f to estimate the DQ approximation of ||f^ l * z || (see the two 
inequalities above equation (29a)). 

Thus, the analysis in f!3f suggests that the estimates for \\r] R ' ltz \\ and ||?7^ z * z || are 
suboptimal with respect to r. The numerical tests in Section^ of this report, however, 
suggest that these estimates are actually optimal, just as in the finite element case. In 
the analysis that follows, we will use the insight from the numerical results in Section^ 
and assume that the Ritz projection yields optimal error estimates with respect to r 
in the L 2 -norm: 

h mtz \\ <C||r?™* erp || , (3.27) 
\\v? itz \\<C\\v interp \\, (3.28) 

We emphasize again, however, that the analysis in fTSj/ does not provide theoretical 
support for (f3~27| - f£2g|) . 

3.2. Case II (V = V no - DQ ). This approach was used in gHH]. The motivation 
for this approach is the following: In Case I (V = V D Q), the first term on the RHS 
of (|3.16|) . {rjt,v r ), yields a term ||?7t|| that stays in all the subsequent inequalities, 
including the final error estimate (|3.24[) . Chapelle et al. proposed in [4] a different 
approach that eliminated the (qt,v r ) term in p,16|) . Their approach was straightfor- 
ward: Instead of using the Ritz projection (as in Case I), they used the L 2 projection. 
That is, they chose w r := P r (u), where P r (u) is the L 2 projection of u, given by 

(u-P r (u),v r ) = Vv r eX r . (3.29) 

Next, we show how the error analysis in Case I changes with w r = P r (u) as in [4]. 
Using ([3T21?]) , (j3~TTj) becomes 

((j> r ,t,v r ) + v(V(t> r ,Vvr) =J^?r+^(Vr?,Vu r ). (3.30) 
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Wc emphasize that it is the cancelation of the first term on the RHS of (|3.30[) that 
yields error estimates that do not require the DQs. We let v r := <f> r in (|3.30|) and we 
apply Cauchy-Schwarz inequality to the remaining term on the RHS: 

\ ^!l<MI 2 + H|v<M| 2 <v||v»?|| l|v<M|. (3.3i) 

The error analysis can then proceed in several directions. 

3.2.1. Approach II. A. One approach is to use Young's inequality in (|3.31l) to 

get 

\ |ll<M| 2 + v l|V0,|| 2 < \ llV/yll 2 + \ ||V«M| 2 , (3.32) 

which implies 

^II^H 2 + ^l|V0 r || 2 <^||Vr ? || 2 . (3.33) 
Noticing that the second term on the LHS of (|3 . 33[) is positive, we get 

|||0r|| 2 <H|Vr/|| 2 . (3.34) 

Comparing the estimate p.34[) with the estimate (|3.20p in Case I, we notice that the 
former has || V77H on the RHS, whereas the latter has \\i]t\\. Thus, using inequalities 
(|3.26p - p.28p in Remark l3.2[ we conclude that Approach II. A will yield error estimates 
that are suboptimal with respect to r. 

We also note in passing that estimate (|3.34p suggests that Approach II. A will 
yield error estimates that are suboptimal with respect to h as well. 

3.2.2. Approach II. B. The other way of continuing from (|3.31|) is to apply the 
POD inverse estimate (|3.5j) : 

||V<M| <C inv (r)\\(j> r \\, (3.35) 
where C inv (r) := y/\\S r \\ 2 . Using (l3~35|) in (|33D yields 

\ j^A 2 + v ||V0 r || 2 < C mv {r) v ||Vi7|| ll^ll- (3.36) 

Dropping v ||V</v|| 2 in (|3 . 36[) and simplifying the resulting inequality by \\<fir\\ (as we 
did in (|3"3Tj) ). we get 

™H0r||<C in „(r)i/||V77||. (3-37) 

Comparing estimate (|3 . 3T[) with estimate ([3734]) in Approach II. A, we note that both 
estimates have ||V?7|| on the RHS. In addition, estimate p.37|) has Ci nv (r) on the 
RHS, which increases the suboptimality with respect to r (see Remark 13. 1| . Thus, 
estimate (|3.37[) suggests that Approach II. B yields estimates that are suboptimal with 
respect to r (and h), just as Approach II. A. 
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3.2.3. Approach II. C. Since both Approach II. A and Approach II. B yield error 
estimates that are suboptimal with respect to r in the L 2 -norm, one can try instead 
to prove optimal error estimates in the i^-norm. To this end, we use the approach 
in |23j and, instead of choosing v r := (f> r in (|3.30[) . we choose v r := <pr.t- 

Ur.t\\ 2 + \ ^H V ^'!| 2 ^ V H Vr ?H HV^r,t||. (3-38) 

Applying Young's inequality and the POD inverse estimate (|3.35|) in (|3.38j) . we get 
ll^|| 2 + ^||V0 r || 2 < H|V»?|| ||V0 r , t || 

< yCH 2 ||Vr;|| 2 + 2 ^ (r)2 l|V0 r , t || 2 



which implies 



t!Q5l ,,2 1 

< T a™(r) 2 ||V7 ? || 2 + -||0 ) , f || 2 , (3.39) 



d ||V0 r || 2 <vC mv (r) 2 \\V V \\ 2 . (3.40) 



dt 



In contrast with estimate (|3.34p in Approach II. A and estimate (|3.37p in Approach 
II. B, estimate (|3.40j) seems to yield error estimates that are optimal with respect to 
r. As in estimates (|3.34[) and p.40p . estimate (|3.40p contains the term ||V?7|| on the 
RHS. We note, however, that this term does not cause any problems, since now we 
are considering the T^-seminorm of the error. The factor Ci nv {r) 2 in p.40p . however, 
increases the suboptimality with respect to r (see Remark l3.ip . Thus, estimate (|3.40l) 
suggests that Approach II. C yields estimates that are suboptimal with respect to r), 
just as Approaches II. A and II. B. 

Since for Case I (V = V D ®) in Section [3~T1 we did not prove error estimates in 
the Tif ^norm, for a fair comparison with Approach II. C, we prove these error estimates 
below. To this end, we let v r := <j) r _t in (|3.16p : 

{<t>r,u4>r,t) + V (V^ r , V<?V,t) = (Vt ^r.t)- (3-41) 

Applying Young's inequality on the RHS of (|3.4ip . we get 



Ml 2 



^l|v0 r || 2 <||^j|||<M| 



which implies 



2 dt' 



-I 
dt 



<\ht\\ 2 + \Ur,.\\\ (3.42) 



V<M| 2 <-|M| 2 . (3.43) 



Comparing (|3.43p with (|3.40p in Approach II. C, we note that the latter does not 
contain the factor Ci nv (r) 2 . Thus, Case I (V = V D ®) in Section |3"7T1 yields optimal 
error estimates with respect to r, as opposed to Approach II. C. 

We also note that, at first glance, estimate (|3.43[) suggests that one can get 
superconvergence in the i/^-norm. As mentioned in [23], however, when applying 
the triangle inequality one obtains the expected convergence rate: 

l|Ve||<||V77|| + ||V&.||. (3-44) 
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3.2.4. Approach II. D. Approaches II. A, II. B, and II. C suggest that the point- 
wise in time (i.e., in the C°(0,T; 7J 1 (0))-norm) error estimates are suboptimal with 
respect to r. Thus, we derive error estimates in the solution norm (i.e., in the 
L 2 (0, T; H x (f2))-norm) . Integrating (j3~33)) from to T, we get 

Ur\\ 2 {T) + V - J T ||V0 r || 2 (.s)ds < ||0 r || 2 (O) + V -j ||Vr ? || 2 ( S ) da. (3.45) 

Estimate in (|3.45[) is optimal with respect to r. 

We note that Proposition 3.3 in [4] yields a similar estimate. 

Case II (i.e., V = v no - D< ^) yields the following general conclusions when the L 2 
projection is used instead of the standard Ritz projection: If the error is computed 
pointwise in time (i.e., in the C°(0, T; L 2 (fi))-norm and the C°(0, T; i/ 1 (f2))-norm), 
then the error estimates are suboptimal with respect to r. This is the main message 
of Approaches II. A, II. B, and II. C. If, however, the error is computed in the solution 
norm (i.e., in the L 2 (0, T; H 1 (f2))-norm), then the error estimates are optimal with 
respect to r. This is the main message of Approach II. D. 



4. Numerical Results. The main goal of this section is to numerically inves- 
tigate the rates of convergence with respect to r of the POD-G-ROM (|2.11j) in the 
two cases considered in Section [3J Case I (V = V DQ ) and Case II (V = V no - DQ ). 
Although the error analysis in Section [3] has been centered around the (linear) heat 
equation, in this section we consider both the heat equation (Section 14. ip and the 
nonlinear Burgers equation (Section 14. 2|) . 

To measure the errors in the two cases (i.e., V — V D< ^ and V — v no - D Q), the 
same norms as those used in Section [3] are used in this section. Denoting the error at 
time tj by e j : = u r h (-,tj) — u(-,tj), the following norms are considered : the error in the 
C°(0, T; L 2 (J7))-norm, approximated by £c°(L 2 ) = m -^y ll e ilU 2 (0); the error in the 

C°(0, T; i7 1 (f2))-norm, approximated by £c"(H 1 ) = Il e j'llffi(f2) j an d the error in 

the L 2 (0,T;i7 1 (^))-norm, approximated by £ L 2 {H i) = Ij^rj X) ll e jllj?ifnv For 

0<j<N 



d 

clarity, we also use the following notation: A r = a X ■ 

As mentioned at the beginning of Section [31 the POD-G-ROM (|2.1ip error esti- 
mates are optimal if the following statements hold: 

(i) The L 2 -norm of the error scales as L 2 -norm of the POD interpolation error. 
Using (|3.ip and (|2.8l) - (|2.9p . this statement is equivalent to 




£co iL i) =0\, \ Xj , (4.1) 



(ii) The TJ^^-norm of the error scales as iJ^^-norm of the POD interpolation error. 
Using ([372]>. ([33]) . and ((2~8l) - ((2~9|) . this statement is equivalent to 



S 



(i 



\Sdh E A H • ( 42 ) 

j=r+l 
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Based on the error analysis in Section [31 we expect the following convergence 
rates with respect to r: 



Table 4.1: Theoretical convergence rates for the no-DQ and the DQ cases. 





no.DQ 


DQ 


£c°(L 2 ) 


suboptimal 
(13.341): (13.371) 


optimal 
(|3.24|) 




suboptimal 
(3.40) 


optimal 
(|3.44|) 




optimal 
(13.45} 


optimal 



4.1. Heat Equation. We consider the one-dimensional heat equation (|2.1[) with 
a known exact solution that represents the propagation in time of a steep front: 



u{x, t) = sin(7ra;) 



1 1 






1 


— arctan 








7T \ 






+ 2 



(4.3) 



where x G [0, 1] and i £ [0, 1]. The constant c in (|4.3|) controls the steepness of the 
front. In all the numerical tests in this section, we used the value c = 100. The value 
of the diffusion coefficient used in the heat equation (|2.1[) is v = 10~ 2 . Piecewise 
linear finite elements are used to generate snapshots for the POD-G-ROM (|2.11l) . A 
mesh size h = 1/1024 and the Crank-Nicolson scheme with a time step At = 10~ 3 
are employed for the spatial and temporal discretization. The time evolution of the 
finite clement solution is shown in Figure |4~T1 In total, 1001 snapshots are collected 
and used for generating POD basis functions in L 2 space. The same numerical solver 
as that used in the finite element approximation is utilized in the POD-G-ROM. 

By varying r, the number of basis functions used in the POD-G-ROM, we check 
the rates of convergence with respect to r for the noJDQ and the DQ cases. The 
errors are listed in Table 14721 (in the no_DQ case) and in Table [4731 (in the DQ case). 
To visualize the rates of convergence with respect to r, these errors with their linear 
regression plots are drawn in Figure 14.21 The convergence rate of the error in the 
C°(L 2 )-norm, £ c o(£2), is supcroptimal in the DQ case and suboptimal in the no-DQ 
case. This supports the theoretical rates of convergence proved in Section [3] and sum- 
marized in Table 14.11 although the suboptimality in the no_DQ case is mild. The 
convergence rate of the error in the C°(7J 1 )-norm, fc°(ff 1 )i is superoptimal in the 
DQ case and optimal in the noJDQ case. The optimality in the no_DQ case seems to 
contradict the theoretical rates of convergence proved in Section [3] and summarized in 
Table |4~T1 The convergence rate of the error in the L 2 (i/ 1 )-norm, £ L 2^ H i^ : is superop- 
timal both in the DQ case and in the no-DQ case. This supports the theoretical rates 
of convergence proved in Section [3] and summarized in Table 14.11 Although we seem 
to get optimal rates of convergence for £qo in the noJDQ case, we emphasize that 
the convergence rates in the DQ case in all three norms are much higher than (and 
almost twice as high as ) the corresponding rates of convergence in the no-DQ case. 
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Fig. 4.1: Heat equation. Fine resolution finite element solution used to generate the 
snapshots. 



Table 4.2: Heat equation. Errors in the noJJQ case. 



r 


A r 


£c°(L 2 ) 


&C a (H 1 ) 




3 


5.72c-02 


9.46e-02 


2.30e+00 


1.59c+00 


5 


2.71c-02 


4.70e-02 


1.58c+00 


1.14e+00 


7 


1.58c-02 


3.69e-02 


1.38c+00 


8.22e-01 


10 


7.34e-03 


1.57c-02 


8.54e-01 


5.31e-01 


13 


3.84e-03 


7.78c-03 


5.84e-01 


3.50e-01 



4.2. Burgers Equation. In this section, we consider the one-dimensional Burg- 
ers equation. As mentioned at the beginning of Section |4j the error estimates proved 
in Section [3] are valid for the (linear) heat equation, but not necessarily valid for the 
nonlinear Burgers equation. Nevertheless, to gain some insight into the range of valid- 
ity of the theoretical development in Section^ we investigate the rates of convergence 
with respect to r in the no_DQ and the DQ cases for the nonlinear Burgers equation: 



u(x, 0) = 
u(x, t) = 



x+UUa 

u {x) 
g(x,t) 



f 



in SI x (0, T] , 
in SI , 

on dSl x (0, T] 



(4.4) 



The initial condition is 



u (x) 



if x € (0, |] 

if x e (|,l) , 



(4.5) 



which is similar to that used in |13j . The diffusion parameter is v = 10 2 , the forcing 
term is / = 0, O = [0,1], and T = 1. The boundary conditions are homogeneous 
Dirichlet, that is, u(0,t) = u(l,t) = for all t £ [0, 1]. 

To generate snapshots, we use piecewise linear finite elements with mesh size 
h = 1/1024 and the backward Euler method with a time step At = 10~ 4 , and save 



16 



TRAIAN ILIESCU AND ZHU WANG 
Table 4.3: Heat equation. Errors in the DQ case. 



r 


A r 


£c°(L 2 ) 




^L 2 {m) 


19 


5.49c-02 


7.15e-03 


4.96c-01 


3.19c-01 


23 


2.95e-02 


2.03e-03 


1.98e-01 


1.19e-01 


28 


1.41e-02 


6.52e-04 


7.97e-02 


4.91e-02 


33 


6.75e-03 


2.41e-04 


3.68e-02 


2.80e-02 


37 


3.76e-03 


8.60e-05 


2.67e-02 


2.29e-02 




Fig. 4.2: Heat equation. Plots of errors in C°(L 2 )-norm (top, left), C°(iJ 1 )-norm 
(top, right), and L 2 (L 2 )-norm (bottom). 



data at each time instance. The time evolution of the finite clement solution is shown 
in Figure 14.31 All snapshots are used for the POD basis generation and the same 
numerical solver is used in the POD-G-ROM. 

By varying r, we check the rates of convergence with respect to r for the no_DQ 
and the DQ cases. Since the exact solution of the Burgers equation (|4.4[) with the 
initial condition (|4.5[) is not known, we consider the errors between the POD-G-ROM 
results and the snapshots. The errors are listed in Table 14.41 (in the no-DQ case) 
and in Table l4~5l (in the DQ case). These errors with their linear regression plots 
are drawn in Figure l4~4l The convergence rate of the error in the C°(L 2 )-norm, 
£c°(l 2 )i is superoptimal in the DQ case and strongly suboptimal in the no_DQ case. 
This clearly supports the theoretical rates of convergence proved in Section [3] and 
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o o t 



Fig. 4.3: Burgers equation. Fine resolution finite element solution used to generate 
the snapshots. 



summarized in Table |4~T1 The convergence rate of the error in the C°(i? 1 )-norm, 
£(70( ff i), is superoptimal in the DQ case and extremely suboptimal in the no-DQ 
case. This strongly supports the theoretical rates of convergence proved in Section [3] 
and summarized in Table |4~T1 Finally, the convergence rate of the error in the L 2 (7J 1 )- 
norm, £j,2/#i), is superoptimal both in the DQ case and in the no-DQ case. This 
supports the theoretical rates of convergence proved in Section [3] and summarized in 
Table 14.11 The convergence rates in all three norms clearly support the theoretical 
rates of convergence proved in Section [3] and summarized in Table 14.11 We also 
emphasize that the convergence rates in the DQ case in all three norms are much 
higher than ( and at least twice as high as ) the corresponding rates of convergence in 
the no-DQ case. 

Table 4.4: Burgers equation. Errors in the no-DQ approach. 



r 


A r 


£c°(L 2 ) 






3 


8.74c-02 


2.38c-01 


4.48c+01 


2.34c+00 


5 


3.95e-02 


1.60c-01 


4.44c+01 


1.76c+00 


7 


1.97c-02 


1.17e-01 


4.37c+01 


1.24c+00 


9 


1.02e-02 


9.05e-02 


4.28e+01 


8.94e-01 


11 


5.47e-03 


7.01e-02 


4.14c+01 


6.84e-01 



5. Conclusions. The effect of using or not the snapshot DQs in the genera- 
tion of the POD basis (the DQ and the no_DQ cases, respectively) was investigated 
theoretically and numerically. The criterion used in this theoretical and numerical 
investigation was the rate of convergence with respect to r of the POD-G-ROM solu- 
tion to the exact solution, where r is the number of POD basis functions used in the 
POD-G-ROM. 

The error estimates in Section [3] yielded the following conclusions: In the DQ 
case, the convergence rates were optimal in all three norms considered (the C°(L 2 )- 
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Table 4.5: Burgers equation. Errors in the DQ approach. 



r 


A r 


£c°(L 2 ) 






16 


1.29c-01 


6.83e-03 


1.18c+00 


4.66c-01 


19 


6.95e-02 


2.99e-03 


4.11e-01 


2.22e-01 


22 


3.69e-02 


1.31e-03 


1.98e-01 


1.08e-01 


25 


1.92e-02 


4.29e-04 


9.74e-02 


5.29e-02 


27 


1.23e-02 


1.88e-04 


6.02e-02 


3.27e-02 



o no-DQ 

- - LR order 0.44 

□ DQ 

— LR order 1.26 




,/||S,||Ej„. + 1 A 



Fig. 4.4: Burgers equation. Plots of errors in C°(L 2 )-norm (top, left), C°(i? 1 )-norm 
(top, right), and L 2 (L 2 )-norm (bottom). 



norm, the C°(i/ 1 )-norm and the L 2 (L 2 )-norm. In the no-DQ case, the convergence 
rates were suboptimal in the C°(L 2 )-norm and in the C°(i/ 1 )-norm, and optimal in 
the L 2 (L 2 )-norm. 

The numerical results in Section |4] for the (linear) heat equation and the (non- 
linear) Burgers equation confirmed the conclusions suggested by the theoretical error 
estimates in Section [3] In the DQ case, the convergence rates were superoptimal in 
the C°(L 2 )-norm, the C°(iJ 1 )-norm, and the L 2 (L 2 )-norm. In the noJJQ case, the 
convergence rates were suboptimal in the C°(L 2 )-norm and in the C°(7J 1 )-norm, and 
optimal in the L 2 (L 2 )-norm. The only departure from the theoretical conclusions was 
that, in the noJJQ case, the convergence rate in the C°(7J 1 )-norm for the heat equa- 
tion was optimal. We note, however, that, in the no-DQ case, the convergence rate 
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in the C> (H )-norm for the Burgers equation was suboptimal, as predicted by the 
theoretical results in Section [3J We also emphasize that, for both the heat equation 
and the Burgers equation, the convergence rates in the DQ case in all three norms 
were much higher than (and usually at least twice as high as) the corresponding rates 
of convergence in the no-DQ case. 

The theoretical error estimates in Section|3]and the numerical results in Scction|4] 
strongly suggest the following conjecture: "The snapshot DQs should be used in the 
generation of the POD basis in order to achieve optimal rates of convergence with 
respect to r, the number of POD basis functions utilized in the POD-G-ROM. We 
also conjecture that using the snapshot DQs in the generation of the POD basis 
could alleviate some of the degrading of convergence with respect to r seen in, e.g., 
PS El [TJ [2] . We intend to investigate this conjecture in a future study. 
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